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Description  of  Problem  and  Approach 

Monte  Carlo  approaches  to  solving  problems  with  many  degrees  of  freedom 
are  a  class  of  statistical  methods  having  in  common  the  generation  of  "random” 
numbers.  In  the  past  few  years,  Monte  Carlo  approaches  have  seen  increased 
application  in  a  number  of  diverse  fields.  In  particular,  quantum  mechanical 
Monte  Carlo  (QMC)  methods1'13  have  been  successfully  used  for  the  treatment  of 
molecular  problems3-58'12.  What  we  mean  here  by  QMC  is  a  Monte  Carlo  pro¬ 
cedure  which  solves  the  Schrodinger  equation  This  is  to  be  distinguished  from 
so-called  variational  Monte  Carlo,  in  which  one  obtains  expectation  values  for  a 
given  trial  wave  function. 

This  ability  to  stochastically  solve  the  Schrodinger  equation  provides  an 
alternative  to  conventional  techniques  of  quantum  chemistry.  Early  work8  has 
shown  that  highly  accurate  total  energies  and  correlation  energies  can  be 
obtained  by  QMC.  In  fact,  in  a  procedurally  simple  manner,  accuracies  exceed¬ 
ing  those  of  the  best  06  initio  configuration  interaction  calculations  have  been 
obtained. 

Much  of  chemistry  takes  place  predominantly  in  the  valence  electrons  of  a 
system.  Thus  the  quantities  of  interest  are  usually  small  differences  of  large 
total  energies.  If  QMC  is  to  be  useful  in  calculating  binding  energies,  affinities, 
reaction  barriers,  etc.,  it  needs  not  only  to  be  able  to  calculate  accurate  total 
energies,  but  also  these  more  relevant  energy  differences .  This  is  a  far  more 
difficult  task  for  Monte  Carlo,  since  a  statistical  uncertainty  of  as  little  as  0.1%  in 
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the  separate  total  energies  can  mask  the  sought-after  energy  difference.  To 
reduce  the  statistical  error  to  the  level  needed  by  "brute  force”  is  costly  in 
computer  time,  as  the  standard  deviation  decreases  only  as  ( CPU  time)' #. 
Algorithmic  developments,  such  as  differential  QMC14,  hold  promise  for  reduc¬ 
tions  in  variance  through  correlated  sampling  techniques.  Another  approach  is 
based  on  noting  that  the  variance  decreases  as  ♦ t  better  approximates  the  true 
eigenfunction.  To  take  advantage  of  this  we  have  developed  an  iterative  pro¬ 
cedure  for  improving  'Iff.  We  describe  this  approach  in  the  next  section. 

Also  under  development  is  the  use  of  pseudo-potentials  in  QMC  to  describe 
the  core  electrons;  in  this  way  only  the  valence  electrons  need  to  be  involved  in 
the  QMC  calculation.  One  wants  also  to  be  able  to  use  QMC  to  calculate  potential 
energy  surfaces,  excited  states,  dipole  moments,  and  other  molecular  proper¬ 
ties.  Research  is  continuing  along  these  lines.  We  describe  in  the  next  section 
some  of  the  results  we  have  obtained.  Among  these  results  is  the  first  calcula¬ 
tion  of  an  excited  state  energy  by  Monte  Carlo.  Before  going  into  more  detail,  in 
what  follows  we  describe  the  actual  QMC  approach. 

Briefly,  the  procedure  is  to  simulate  the  quantum  system  by  allowing  it  to 
evolve  under  the  time-dependent  Schrodinger  equation  in  imaginary  time.  It  is 
easy  to  show0  that  the  use  of  imaginary  time  causes  the  system  to  approach  a 
stationary  state  which  is  the  lowest  state  of  a  given  symmetry.  Properties  may 
then  be  "measured"  as  averages  over  the  resulting  equilibrium  distribution. 

By  writing  the  imaginary-time  Schrodinger  equation  with  a  shift  in  the  zero 
of  energy  as 

+  [ET-V(E))*(R.t)  ,  (1) 

we  see  that  it  may  be  interpreted  as  a  generalized  diffusion  equation.  The  first 
term  on  the  right-hand-side  is  the  ordinary  diffusion  term,  while  the  second 
term  is  a  position-dependent  rate  (or  branching)  term.  For  an  electronic  sys¬ 
tem,  D=HZ/  2m,,  ^  is  the  three-N  dimensional  coordinate  vector  of  the  N  elec¬ 
trons,  and  V(fi)  is  the  Coulomb  potential.  Since  diffusion  is  the  continuum  limit 


of  a  random  walk,  one  may  simulate  Eq.  (1)  with  the  function  'f'  (note,  not  4'8)  as 
the  density  of  “walks".  The  walks  undergo  an  exponential  birth  and  death  as 
given  by  the  rate  term.  This  connection  between  a  quantum  system  and  a  ran¬ 
dom  walk  was  first  noted  by  Metropolis,  who  attributes  it  to  Fermi15. 

The  steady-state  solution  to  Eq.  (1)  is  the  time-independent  Schrodinger 
equation.  Thus  we  have  ^(R.tJ^tfifR),  where  ^  is  an  energy  eigenstate.  The 
value  of  Et  at  which  the  population  of  walkers  is  asymptotically  constant  gives 
the  energy  eigenvalue.  Early  calculations  employing  Eq.  (1)  in  this  way  were 
done  by  Anderson  on  a  number  of  one-  to  four-electron  systems3. 

Unfortunately,  in  order  to  treat  systems  larger  than  two  electrons,  the 
Fermi  nature  of  the  electrons  must  be  taken  into  account.  The  antisymmetry  of 
the  eigenfunction  implies  that  '{'  must  change  sign;  however,  a  density  (e.g.  of 
walkers)  cannot  be  negative.  To  handle  this,  Anderson  made  simplifying 
assumptions  about  the  positions  of  the  nodes.  His  method  is  ad  hoc ,  and  not 
readily  generalizable.  Another  method  which  imposes  the  antisymmetry,  and  at 
the  same  time  provides  more  efficient  sampling  (thereby  reducing  the  statisti¬ 
cal  "noise"),  is  importance  sampling  with  an  antisymmetric  trial  function  '{'r 
(see  e.g.  Ref.  B).  The  zeroes  (nodes)  of  become  absorbing  boundaries  for  the 
diffusion  process,  which  maintains  the  antisymmetry.  The  additional  boundary 
condition  that  4'  vanish  at  the  nodes  of  '{'r  is  the  fixed-node  approximation.  The 
magnitude  of  the  error  thus  introduced  depends  on  the  quality  of  the  nodes  of 
+ t(R).  and  vanishes  as  ♦p  approaches  the  true  eigenfunction.  To  the  extent 
that  'Vf  is  a  good  approximation  of  the  wave  function,  the  true  eigenfunction  is 
almost  certainly  quite  small  near  the  nodes  of  Thus  one  expects  the  fixed- 
node  error  to  be  small  for  reasonable  choices  of  'frp.  Work  on  a  number  of  sys¬ 
tems  has  borne  this  out8,9,11.  In  addition,  this  error  is  variationally  bounded. 

To  implement  importance  sampling  and  the  fixed-node  approximation,  Eq. 
(1)  is  multiplied  on  both  sides  by  'frp,  and  rewritten  in  terms  of  the  new  probabil¬ 
ity  density  /  The  resultant  equation  for  f(R,t)  may  be 


written 


f£-=  DV*f  +[ET-EL(&)]f  -DHfF<l(&)]  ■  (2) 

The  local  energy  El(R),  and  the  "quantum  force"  Fq{R)  are  simple  functions  of 
'fry  given  by 

£xCff)=//*7't2)/*rCff)  .  (3a) 

and 

/Vtff)S2V*7.Cff)/*rCfi)  .  (3b) 

Equation  (2),  like  Eq.  (1)  is  a  generalized  diffusion  equation,  though  now  with  the 

addition  of  a  drift  term  due  to  the  presence  of  Fq. 

In  order  to  perform  the  random  walk  implied  by  Eq.  (2)  we  use  a  short-time 
approximation  to  the  Green's  function  which  is  used  to  evolve 
f  (R.t +t).  This  process  is  iterated  to  large  t.  Such  an  approach 
becomes  exact  in  the  limit  of  vanishing  time-step  size,  r. 

Progress  During  Current  Year 

(1)  We  have  developed  a  method  of  improving  an  importance  function  through  an 
iterated  use  of  wave-function  scaling.  For  use  in  QMC  one  wants  a  trial  function 
which  is  as  simple  as  possible,  since  it  will  require  repeated  evaluation  at  each 
step  of  the  random  walk.  Yet  one  wants  a  function  which  provides  accurate 
results.  In  principle,  since  QMC  solves  the  Schrodinger  equation,  one  should 
obtain  accurate  results  regardless  of  the  choice  of  if?-  However,  as  we  noted 
already,  for  Fermi  systems  inaccurate  nodes  in  ♦ t  will  lead  to  a  small  error 
when  the  fixed-node  approximation  is  used.  Furthermore,  the  statistical 
"noise”  will  be  large  for  a  poor  choice  of  <'f. 

We  have  found  in  our  work8,11  that  a  single  determinant  'f'r  with  only  a 
double-zeta  basis  set  places  the  nodes  extremely  well  as  determined  by  the 
quality  of  the  computed  total  energies.  Increasing  the  basis  set  beyond  double 
zeta  appears  to  offer  insignificant  gain  in  either  accuracy  (i.e.  the  fixed-node 
error  does  not  noticeably  decrease)  or  precision  (the  statistical  uncertainty,  for 


equal  computing  time,  remains  essentially  unchanged).  In  practice  we  have 
included  an  electron-electron  Jastrow  factor  in  our  functions  'f'r  in  order  to 
reduce  statistical  fluctuations,  and  in  some  cases  we  have  also  included  an 
electron-nuclear  factor.  Neither  factor  affects  the  positioning  of  the  nodes,  and 
hence  the  fixed-node  energies. 

Since  the  fixed-node  QMC  approach  involves  an  approximation  in  the  place¬ 
ment  of  the  nodes,  and,  in  addition,  in  many  applications  the  statistical  uncer¬ 
tainty  needs  to  be  further  reduced,  we  have  developed  an  iterative  approach  for 
globally  (rather  than  locally)  performing  this  correction. 

Note  first  that  the  additional  boundary  condition  that  the  eigenfunction 
vanish  at  the  nodes  of  ♦ r  will  generally  give  a  solution  which  fails  to  satisfy  the 
virial  theorem.  Thus  the  fixed-node  expectation  value  of  V  is  not  exactly  —2 T. 
(Here  we  assume  an  equilibrium  geometry.)  We  may  therefore  consider  the 
fixed-node  eigenfunction,  which  we  shall  denote  ^(foj),  as  a  variational  function 
which  may  be  further  optimized  by  scaling16.  In  our  notation,  the  caret  indi¬ 
cates  that  <p  carries  the  fixed-node  constraint,  and  |rt{  indicates  the  set  of  all 
coordinates. 


Let  us  define 

V(l)  ■  <£0rJ)|  V|£(|r*})>  . 

(4a) 

r(i)  =  <?«rt0in?(ta»>  . 

(4b) 

and  the  scaled  quantities  1^(77)  and  T(rj)  analogously  in  terms  of  the  scaled  func¬ 
tion  The  expressions  in  Eq.  (4)  must  of  course  be  divided  by 

<£(lr«{)l£Or<i)>  if  9  is  not  normalized.  It  is  readily  established  that 
V(r))=T)V(i)  since  the  Coulomb  potential  scales  as  1/r.  Similarly,  T(r))=r)zT(l) 
since  V2  scales  as  1/r2.  Combining  these  expressions  one  obtains 

E(r,)=r)V(l)+ri*T(l)  (5) 

Varying  Eq.  (5)  with  respect  to  77  minimizes  E(r})  at 

77=-K(l)/2r(l) 


and 


(6a) 


£(*))="  HO2/ 47(1)  (6b) 

Thus  the  function  £(»j$rt  j)  has  a  lower  variational  energy  than  and  in 

addition  satisfies  the  virial  theorem  since  —V(r))/ ZT(r))  =  -r)~lV(l)/  2T(l)  =  1. 
Note  that  the  global  scaling  has  uniformly  expanded  or  contracted  the  nodal 
surfaces  originally  present  in  if.  As  we  demonstrate  below,  these  new  nodes  are 
better  than  the  original  nodes  of  'if.  However,  is  no  longer  an  eigen¬ 

function  of  the  Schrodinger  equation.  Thus  we  may  iterate  the  above  procedure 
starting  with  the  new  nodes—i.e.  using17  a  4'j'  whose  nodes  are  those  of  ^(Tjfo}) 
(see  Fig.  1.).  Such  a  function,  ift\  may  be  obtained  by  replacing  all  coordinates 

in  if  by  (Essentially  this  involves  scaling  all  the  orbital  exponents 

and  the  inter-atomic  separations.)  Now  starting  with  the  QMC  method  con¬ 
verges  to  an  eigenstate  ^,*(T?lri j).  Because  ^lKr}[ril)  has  the  same  nodes  as 
?(?)lrtO«  must  have  the  lower  energy  since  it  is  the  exact  solution  for  these 
nodes.  Again,  due  to  the  fixed-node  boundary  condition  with  the  new  nodes,  the 
virial  theorem  may  not  be  satisfied,  resulting  in  rj'  =  -V/2.T  / 1  (for  Thus 

we  rescale  by  r\’  to  obtain  which  has  a  lower  variational  energy  and 

again  satisfies  the  virial  theorem.  The  expanded  or  contracted  nodes  may  then 
be  fed  back  into  a  ift  and  the  process  repeated.  It  is  expected  that  the 
sequence  rj,  if,  r}",  ...  rapidly  converges  to  unity,  so  that  no  appreciable  gains 
will  be  obtained  beyond  the  first  few  iterations.  Figure  1  gives  a  schematic  illus¬ 
tration  of  this  iterative  procedure.  Since  the  fixed-node  energies  for  the 
sequence  of  functions  ^(foj),  SP^IHTJtr*  1),  ?^(T)’T)lrtD—  is  of  decreasing  energy, 
the  nodes  improve  upon  scaling. 

(2)  We  have  calculated  the  energy  of  the  first  excited  state  of  methylene18  in 
order  to  obtain  the  (until  recently)  elusive  singlet-triplet  splitting.  This  is  the 
first  QMC  calculation  of  an  excited  state.  Our  results  are  in  excellent  agreement 
with  the  most  recent  experiments. 

Also,  as  shown  in  Table  I,  the  QMC  total  energies  for  both  the  singlet  and  the 
triplet  states  of  methylene  compare  favorably  with  Cl  calculations.  For  the  best 
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trial  function  used,  the  total  energy  is  correct  to  better  than  0.008  h  (5 
kcal/mole)  of  experiment,  or  to  1  part  in  5000.  The  statistical  uncertainty  is 
roughly  half  this  value  (2-2.5  kcal/mole).  The  remaining  error  may  be  attri¬ 
buted  to  the  fixed-node  and  the  short-time  approximations.  This  translates  to  a 
Monte  Carlo  accuracy  of  99.98%  of  the  total  energy  and  96-98%  of  the  correlation 
energy.  Therefore  in  the  present  application,  where  the  time-step  error  is  negli¬ 
gible,  the  fixed-node  error  is  seen  to  be  manageably  small.  Furthermore,  to 
within  1  kcal/mole,  this  error  is  the  same  for  the  two  states.  This  means  not 
only  that  the  absolute  error  is  small,  but  that  there  is  also  a  large  degree  of  can¬ 
cellation  of  this  error  in  evaluating  the  energy  gap.  In  fact,  for  the  energy  gap 
the  error  is  considerably  less  than  the  statistical  uncertainty. 

To  obtain  our  best  estimate  for  T,t  we  calculated  a  weighted  average  of  the 
energy  differences  for  our  various  trial  functions.  Our  final  result  is  TB  =9.4±2.2 
kcal/mole.  This  result  is  in  excellent  agreement  with  the  recent  experimental 
results  of  McKellar  et  al 24 . 

(3)  We  have  calculated  points  along  the  reaction  path  of  the  H  +  Ylg  exchange 
reaction.  Particular  emphasis  has  been  placed  on  the  saddle  point,  for  which 
Liu23  has  performed  the  most  extensive  Cl  calculation  to  date.  Nevertheless,  the 
bound  for  the  barrier  height  which  we  obtained26  by  QMC  is  0.16  kcal/mole  below 
Liu’s  bound  (see  Table  II),  and  probably  lies  within  0.1  kcal/mole  of  the  exact 
answer.  In  addition,  we  were  able  to  obtain  these  results  with  only  single- 
determinant  trial  functions,  and  a  basis  set  expansion  at  only  the  double-zeta 
level. 

The  nodes,  which  are  so  important  in  determining  the  correct  energy,  prove 
to  be  quite  insensitive  to  basis  set  after  the  double  zeta  level.  A  single-zeta  basis 
set,  however,  gives  a  very  poor  nodal  description  (see  Figs.  2-4). 


(4)  Monte  Carlo,  as  alluded  to  earlier,  is  such  a  computationally  intensive 
activity  that  new  techniques  are  needed  if  one  wishes  to  attack  large  systems  or 


obtain  very  high  precision  (e.g.  better  than  99.99%).  One  avenue  we  have 
explored,  in  collaboration  with  the  Advanced  Computer  Architecture  Laboratory 
at  LBL,  is  the  use  of  parallel  computing  architectures.  Briefly,  we  have  found 
that  Monte  Carlo  can  be  readily  made  to  run  at  95%  efficiency  on  an  8  processor 
system.  With  sufficient  memory,  precision  will  scale  as  the  square  root  of  the 
number  of  processors,  while  computing  time  for  a  fixed  precision,  will  scale 
inversely  almost  linearly  with  processors.  We  explored  several  different  direc¬ 
tions  for  parallelizing  QMC,  as  well  as  load-balancing  techniques  to  keep  the 
efficiency  near  100%. 

(5)  An  accurate  calculation  of  the  binding  energy  of  Ng  has  been  a  classically 
difficult  problem  using  traditional  oft  initio  quantum  chemical  approaches. 
Since  the  quantity  of  interest.  Efrmdinf  is  desired  to  better  than  5  kcal/mole  out 
of  a  total  energy  (for  Ng)  of  over  68,000  kcal/mole,  this  is  an  example  where 
QMC  requires  very  high  precision.  Thus  we  have  performed  this  calculation  on 
the  experimental  parallel  processing  system  mentioned  above.  We  have 
obtained  the  Ng  binding  energy  to  be  233  ±  5  kcal/mole.  To  within  the  statisti¬ 
cal  uncertainty,  this  answer  agrees  with  experiment. 

(6)  We  have  also  begun  investigating  the  electron  affinity  of  F,  since  affinities  are 
generally  difficult  by  traditional  approaches.  In  addition,  we  are  exploring  some 
new  directions.  As  mentioned  in  the  previous  section,  molecular  pseudopoten¬ 
tials  hold  promise  for  simplifying  Monte  Carlo  calculations.  Little  loss  of  accu¬ 
racy  is  anticipated.  However,  this  must  be  explored.  More  importantly,  one 
must  learn  how  to  translate  angular  momentum  dependent  “effective-core- 
potentials"  into  a  form  applicable  to  QMC.  Another  direction  we  are  following  is 
the  explicit  calculation  of  derivatives  by  Monte  Carlo.  This  would  enable  direct 
calculation  of  forces  and  more  accurate  determination  of  potential-energy  sur¬ 
faces.  Also,  this  would  provide  the  ability  to  And  equilibrium  geometries. 
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Table  /.  Comparison  of  QMC  results  on  methylene  with  SCF,  Cl,  and  experimen¬ 
tal  values.  The  results  indicated  as  "expt"  are  corrected  for  zero-point  motion 
and  relativistic  effects  to  make  the  comparison  direct.  The  lowest  variance  QMC 
result  for  7*.  is  not  the  difference  of  the  QMC  values  given  for  the  two  states. 
Instead  T,  is  obtained  by  averaging  the  results  of  all  the  'fry's  we  used  in  this 
study. 

Table  II.  Comparison  of  computed  reaction  barrier  heights  for  the  H  +  Hg 
exchange  reaction.  Fixed-node  quantum  Monte  Carlo  provides  an  excellent 
bound,  which  is  better  than  that  given  by  Cl,  and  is  within  0.1  kcal/mole  of  the 


Table  1. 


Method* 

Energy  (hartrees) 
_ 3!i _ _ 

T. 

(Kcal/mole) 

SCF 

-38.9348“ 

-38.S9446 

25.4 

2C-SCF 

— 

-38.9177* 

10.7* 

-39.1160* 

-39.1003* 

9.9* 

CI-SD(Q) 

-39.122° 

-39.105° 

10.7° 

QMC 

-39.140(4) 

-39.128(3) 

9. 4(2.2) 

•■expt" 

-39. 148^ 

-39.133* 

9.55* 

Glossary  of  Methods 

SCF  =  self-consistent  field 
2C-SCF  =  two-configuration  SCF 

2R-C1-SD  =  two-reference  configuration,  single  and  double  excitations  Cl  for 
the  singlet;  one  reference  configuration  CI-SD  for  the  triplet. 

Cl-SD(Q)  =  singles  and  doubles  Cl,  quadruples  estimated. 

QMC  =  quantum  Monte  Carlo  (this  work). 


°  Ref.  19,20. 

6  Ref.  21. 
e  Ref.  19. 

*  using  1C-SCF  for  the  state. 

*  Ref.  22. 

*  Ref.  23. 

9  obtained  by  subtracting  T,  from  the  "expt"  energy  of  B ^  CHg. 
obtained  from  T0  of  McKellar  et.  ai.  (Ref.  24),  corrected  for  zero-point  motion 


and  relativistic  effects. 
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Figure  Captions 

Figure  1.  Schematic  illustration  of  process  for  globally  optimizing  ♦ t •  The  origi¬ 
nal  trial  function  is  'FrO*'*}).  while  the  subsequent  trial  functions  are  4^, 
etc.  The  functions  in  the  middle  column  are  the  solutions  of  the  fixed-node 
Schrodinger  equation,  with  the  nodes  of  the  4'r  to  the  left.  The  final  functions  on 
the  right  are  the  scaled  fixed-node  functions,  but  these  are  no  longer  solutions 
to  the  Schrodinger  equation.  Along  the  path  indicated  by  arrows,  each  function 
has  a  lower  variational  energy  than  that  preceding  it. 

Figure  2.  Exchange  nodes  of  a  single-zeta  quality  trial  function.  The  circles 
represent  the  hydrogen  nuclei.  The  curves  are  cross  sections  through  a  selec¬ 
tion  of  nodal  surfaces  arising  from  the  exchange  antisymmetry.  Full  nodal  sur¬ 
faces  may  be  obtained  by  rotating  the  curves  around  the  internuclear  axis. 
Each  surface  is  obtained  by  fixing  the  position  of  one  electron  on  it  and  finding 
the  locus  of  points  for  the  other  like-spin  electron  at  which  the  trial  function 
vanishes.  It  can  easily  be  shown  that  the  trial  function  is  zero  whenever  both 
like-spin  electrons  are  anywhere  on  this  surface.  Distances  are  in  bohr. 

Figure  3.  Exchange  nodes  for  a  double-zeta  quality  trial  function.  See  Fig.  2  for 
further  explanation  of  this  figure.  Note  how  different  these  nodes  are  from  those 
depicted  in  Fig.  2. 

Figure  4.  Exchange  nodes  for  a  near  Hartree-Fock  quality  trial  function.  See 
Fig.  2  for  further  explanation  of  this  figure.  Note  that  approaching  the  basis-set 
limit  has  little  effect  on  the  nodes,  past  the  single-zeta  level. 
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